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Abstract. We calculate the zero-temperature self-energy to fourth-order perturbation theory in the Hub- 
bard interaction U for the half-filled Hubbard model in infinite dimensions. For the Bethe lattice with bare 
bandwidth W, we compare our perturbative results for the self-energy, the single-particle density of states, 
and the momentum distribution to those from approximate analytical and numerical studies of the model. 
Results for the density of states from perturbation theory at U/W = 0.4 agree very well with those from 
the Dynamical Mean-Field Theory treated with the Fixed-Energy Exact Diagonalization and with the 
Dynamical Density-Matrix Renormalization Group. In contrast, our results reveal the limited resolution of 
the Numerical Renormalization Group approach in treating the Hubbard bands. The momentum distribu- 
tions from all approximate studies of the model are very similar in the regime where perturbation theory 
is applicable, U /W < 0.6. Iterated Perturbation Theory overestimates the quasiparticle weight above such 
moderate interaction strengths. 

PACS. 71.10.Fd Lattice fermion models (Hubbard model, etc.) - 71.27.+a Strongly correlated electron 
systems; heavy fermions - 71.30+h Metal-insulator transitions and other electronic transitions 



1 Introduction 



The Hubbard model serves as a paradigm for strongly cor- 
related electron systems because it combines the two es- 
sential aspects of electrons in solids in a simplistic way. It 
describes spin-1/2 electrons moving on a lattice (band- 
width W) which interact only locally with strength U 
(Hubbard interaction). Despite its structural simplicity it 
is thought to encompass a rich phase diagram. 

At half band-filling, when there is on average one elec- 
tron per lattice site, the Hubbard model contains a zero- 
temperature phase transition from a metal to an insula- 
tor, irrespective of any symmetry breaking jll'2l.'{l lj . In 
the limit of large lattice dimensions the precise na- 
ture of this Mott-Hubbard transition is not entirely clear. 
An exact solution of the problem is not possible, and var- 
ious approximate treatments have led to two conflicting 
scenarios; for a review, see Refs. |2I3| . For more recent 
treatments, see Refs. ()I1 

Discontinuous Transition. 

The gap jumps to a finite value when the density of 
states at the Fermi energy becomes zero at some crit- 
ical interaction strength J7 C| a; the gap is preformed 
above C/ Cj i < U C: 2, and the co-existing insulating state 
is higher in energy than the metal. 



Continuous Transition. 

The gap opens continuously when the density of states 
at the Fermi energy becomes zero, C/ C ,i = U c p = U c . 

This situation calls for the development and appli- 
cation of systematic and controlled techniques such as 
high-order perturbation theory in strong and weak cou- 
pling. Recently, some of us ^2] investigated the Mott- 
Hubbard insulator on a Bethe lattice in the limit of large 
coordination number analytically and numerically. In the 
present work, after the introduction of some definitions 
in Sect. [5J we calculate the one-particle Green function 
of the metallic state at half band-filling up to and in- 
cluding ©[(C//VF) 4 ] in Sect. 21 corrections are of the or- 
der (U/W) 6 due to particle-hole symmetry. Perturbation 
theory to fourth order is found to converge very well at 
U = OAW, but it begins to fail at U » 0.64TF. In Sect.H 
we compare our perturbative results for the density of 
states with those of the Dynamical Mean-Field Theory 
(DMFT), analyzed within two recently developed numeri- 
cal schemes, the Fixed-Energy Exact Diagonalization (FE- 
ED) ^21 and the Dynamical Density-Matrix Renormal- 
ization Group (DDMRG) [HITS] , The DMFT, which be- 
comes exact in the limit of infinite dimensions, requires 
the self-consistent solution of a single-impurity Anderson 
model on n s — ► oo bath sites. We employ the FE-ED for 
n s < 15 and the DDMRG on up to n s = 64 sites. We 
find very good agreement with our perturbative results 
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for U = 0.4W, and can attribute deviations at U = 0.6W 
to the limited accuracy of fourth-order perturbation the- 
ory. Therefore, we conclude that the DDMRG provides a 
reliable description of the correlated metal at all frequen- 
cies. The comparison is less favorable with results from the 
Numerical Renormalization Group (NRG) ^Uj an d Iter- 
ated Perturbation Theory (IPT) which show noticeable 
deviations for intermediate energies, e.g., in the formation 
of the Hubbard bands. 

In Sect.0 we compare our results for the momentum 
distribution n(e) and the quasi-particle weight Z(U) with 
those of the Random Dispersion Approximation (RDA), 
the NRG and IPT. Within the region of validity of our 
perturbation expansion, U < 0.6, the results are almost 
identical except for the IPT which already overestimates 
the quasi-particle weight at weak to moderate interaction 
strengths. Conclusions, Sect.|Bl and two appendices close 
our presentation. 

2 Definitions and basic properties 

In this section, we discuss the basic properties of lattice 
electrons in the limit of infinite dimensions. We define the 
Hubbard Hamiltonian, the one-particle Green function, 
and some related one-particle quantities. 

2.1 Hamilton operator 

We investigate spin- 1/2 electrons on a lattice. Their mo- 
tion is described by 

iJ;o- 

where cf ' , c\. a are creation and annihilation operators for 
electrons with spin a =|,| on site i. Here iij are the 
electron transfer amplitudes between sites i and j, and 
t\\ = 0. Since we are ultimately interested in the Mott- 
Hubbard transition we consider a half-filled band where 
the number of electrons N equals the number of lattice 
sites L exclusively. 

For lattices with translational symmetry, we have ti j = 
t(i— j), and the operator for the kinetic energy is diagonal 
in momentum space, 

T = Y, e ( k ) 6 k,<A* » 

k;<r 

e(k) = i^t(i-j)e- i ( I -j) k . (2) 

The density of states for non-interacting electrons is then 
given by 

P(e) = ^«5(e-e(k)). (3) 

k 

The ra-th moment of the density of states is defined by 

/oo 
dee>(e), (4) 
-CO 



and e = t(0) = 0. 

In the limit of high lattice dimensions and for trans- 
lationally invariant systems, the Hubbard model is char- 
acterized by p(e) alone, i.e., higher-order correlation func- 
tions in momentum space factorize, e.g., |16| 

ftu,q 2 (ei>e2) = -^E5(ei-e(k + q 1 ))(5(e 2 -e(k + q 2 )) 

k 

= P(ei)[<*qi,q!»tf(ei - £ 2) + (1 - V.qsVfe)] . 

(5) 

This observation is the basis for the Random Dispersion 
Approximation (RDA) which becomes exact in infinite di- 
mensions for paramagnetic systems where nesting is ig- 
nored PJE]; see Sect. 03 

In the following, we assume a symmetric bare density 
of states, p{— e) = p(e), of width W. For our explicit cal- 
culations we shall later use the semi-circular density of 
states 

p^)-^fA^f , (H<W72), (6) 

where W = At is the bandwidth. In the following, we shall 
set t = 1 as our energy unit if not otherwise explicitly 
stated. This density of states is realized for non-interacting 
tight-binding electrons on a Bethe lattice of connectivity 
Z — > oo rj7], i.e., each site is connected to Z neighbors 
without generating closed loops, and the electron transfer 
is restricted to nearest-neighbors, = —X/yZ for i and 
j being nearest neighbors and zero otherwise. The limit 
Z — > oo is implicitly understood henceforth. 

The electrons are assumed to interact only locally, and 
the Hubbard interaction reads 




where h^ a = c^ a Ci tCr is the local density operator at site i 
for spin a. This leads to the Hubbard model |18| . 

H = f + UD. (8) 

The Hamiltonian exhibits explicit particle-hole symmetry, 
i.e., H is invariant under the transformation 

PH: ~ (-l)'ci,* ; h,„ -> (-1)^ ■ (9) 

For the Bethe lattice (— l) 1 = +1 for the A sites which 
are surrounded by B sites only (and vice versa) , for which 
(— l) 1 = — 1. A chemical potential p — then guarantees 
a half- filled band for all temperatures 

2.2 Green functions and self-energy 

The time-dependent single-particle Green function at zero 
temperature is given by |19j 

G <T (i,j;t) = -i(f[c i)CT (t)c+J) . (10) 
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Here T is the time-ordering operator, (. . .) implies the With the help of a full set of eigenstates \& n ), we can 

expectation value in the ground state with energy Eo, and write the local Green function in the Lehmann represen- 

(h = 1) tation [Hj 

Ci,a(t) = exp(iHt)ci t<T exp(-iHt) (11) x D (u 1 ) 

is the annihilation operator in the Heisenberg picture. For G a (uj) — J ^du — — — — j^g^^ (^0) 

translationally invariant systems, the Fourier transform of ^ 

the Green function is given by D a {uj > 0) = ^ ~[ X! I (^olci^l^) | 2 S(u> + E Q - E n ) . 

n i 

/°° 1 
Yj e ' ■''C (T (i, j; t) . (12) Consequently, the density of states is obtained from the 
' j imaginary part of the local Green function (|17Jl via 

The Green function in momentum space can be expressed 1 

in terms of the self-energy E a (\t,ui) [32] (17 = 0+), D ^ UJ > = --sgn(w)3G CT (a;) = D a (-u>) . (21) 

1 The latter equality holds due to particle-hole symmetry. 

G CT (k,w) = - - ^ - —^ - i^ S g n ( w ) ' ( 13 ) For a Fermi liquid, JHJ and in (HJ lead to 

In the limit of infinite dimensions, the self-energy is inde- 

£> ff (0)=p(0), (22) 

pendent of momentum 51, -EUfk, to) = E„{u>). As we shall ,, , .. r , „ . . , , . , 

c . r, , ri ir 1 i i i i.e., the density of states at cj = is pinned to its value at 

further discuss m Sect. |3J the self-energy can be calculated [/ — 2"T1 

in a power series in U. As shown by Luttinger 20 , it has , , r , •, , r , JC j 

„.•' 1 he moments M„ of the density of states arc defined 

the following properties at small frequency, 

as 

roc 



, r ,1\ M n = 2 du^D^u). (23) 

uE^to) = [ 1 - — Iw ; 0<cj<cj c , (14) J 

2 n ^ ^ In particular, from l|2(J|) . I|21|) . and the definition of the 

Green function (|10fl one can show that [TiJ\ 



%E a (to) = ; < to < lu c . (15) 



Here < Z < 1, 7 > 0, and a low-energy cut-off lo c <C VF j / 



characterize a Fermi liquid. Mi = — — IEq + U j ■ (24) 

The local Green function G a (u>) is the momentum- ^ ' 

space average of G CT (k, ui), yy e later employ this useful sum rule for calculating 

the ground-state energy. 

G a (uj) = — / dt e iwt G CT (i, i; f) The one-particle spectral function A CT (e; w) is symmet- 

L ~ J -00 ric in to. It is defined by {to > 0) 

= lE^-) (16) A CT (e; W ) = -W -) (25) 

k 7T \ w — e — E a (io) + iT]J 

de P (e) I - 1 



to-e- E a {io) + irysgn(w) ' n ( w _ e _ 9fJ£ CT ( w )) 2 + (9f27 ff (o;) - r?) 2 ' 

Due to particle-hole symmetry, the local Green function Note that 

obeys G CT (i) = — G CT (— t) so that G a {uo) = —G a (—uj), and oo 

thus Ejfjj) = -E a {-w). D a (co) = dep(e)A a {e;to) . (26) 



Because the self-energy only depends on frequency, the 

local Green function can be easily recovered from the self- , , . ,. „ , T , 

a r rm As follows from the Lehmann representation 19 , the spec- 
energy. As seen from HSJ, , , ... . ■ ,!,„ A ,„ , V*' n _ K 



G a {u) = G° (w - E a {u)) . (17) 



tral function is positive semi-definite, A a (e; u>) > 0, so that 
%E a (Lu > 0) < . (27) 



This relation is particularly simple for the semi-circular xhe spectral f unct i n contains a quasi-particle contribu- 
density of states (J6J) where y on r we ight Z(U) near the Fermi energy and an inco- 

herent background contribution, 

( 18 ) A a (e; lj^0) = Z(U)S(u;-Z(U)e)+A^ c (e; to -> 0) . (28) 



with z = o>+ir/sgn(a;). For the Bethe lattice it then follows 
that 



The momentum distribution 







u-E„{u)=G (T {u) + G-\u). (19) n„{e)= \ duA^u) (29) 
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depends on momentum only implicitly via e = e(k). In the 
metallic phase, n a (e) displays a jump discontinuity at the 
Fermi energy, 

n a (e = 0-)-n a (e = 0+) = Z(U) , (30) 

as follows directly from ifTIj). JT5J), (|2*3|) and the fact that 
the self-energy does not depend on momentum. In the 
vicinity of e = 0, the momentum distribution takes the 
form (E = Z(U)e/uj c ) 

„„ (£<<1) = I^ + ^±£<^ ElllE + 0(E) . 

2 ir 

(31) 

This is shown in appendix 1X1 

3 Diagrammatic perturbation theory 

In this section we derive and calculate the diagrams to 
second and fourth order perturbation theory in U. 

3.1 Second order 

The particle-hole transformation J§J can be restricted to 
one spin species only. The Hamiltonian then maps onto 
itself apart from a change in the sign of U. Therefore, the 
Green function obeys 

G„(u;U) = G a {cj;-U) , (32) 

and, correspondingly, the self-energy fulfills 

E„(w;U)=E tr (u;-U). (33) 

Consequently, there are no odd orders in the perturbation 
expansion of the self-energy in U. 

Particle-hole symmetry also guarantees that there are 
no (renormalized) Hartree bubbles. A chemical potential 
1-1 = results in [2] 

lD^> = lI>>o = ~, (34) 

i i 

i.e., the bare Hartree diagrams are exactly canceled by def- 
inition of the Hubbard interaction in JJJ), and, moreover, 
the renormalized Hartree diagram vanishes to all orders 
in perturbation theory. Lastly, there are no Fock contri- 
butions because the Hubbard interaction acts between dif- 
ferent spin species only. 

With these simplifications, only one diagram remains 
in second-order perturbation theory, 



UJ-2 + OJl 




Three independent Green function lines connect the two 
lattice points. Thus, in the limit of infinite dimensions, 
they can be identified with each other and each line 
thus represents a local bare Green function G°(uJi) = 
G a _ a (uii) due to spin symmetry. Note, however, that en- 
ergy conservation must still be obeyed at each vertex. Fol- 
lowing the Feynman rules, the second-order diagram gives 
the contribution |21I22| 

s^liu) = (-i)(i) 2 t/ 2 r ^G°_ a (—i)n>i) (36) 

with the bare polarization bubble 

n» = - f H Gl (a* + u) = ZZS(- W ) . 

(37) 

With the help of the spectral representation (|16(l of the 
local bare Green function 

Gl(u) = / dep(e) — + — , 

JO \<jj-e + 17] UJ + 6 -17] J 

(38) 

a contour integration results in 

fW/2 rW/2 

n° a {uj) = - dei / de 2P (e 1 )p(e 2 ) (39) 

Jo Jo 

( l — !_^) . 

\lj — ei — e 2 + if] ui + ei + e 2 - \7] ) 
Taking the imaginary part gives 

-377° (w > 0) = [ U dep(e)p(cj-e) . (40) 
t Jo 

This representation explicitly shows that the imaginary 
part of the bare polarization bubble vanishes for \u>\ > W . 
This is a consequence of the fact that the bare polarization 
bubble is made up of two bare Green function lines. 

Using the spectral representation of the bare polar- 
ization bubble (|39|l and of the local Green function (|38|l 
in l|36l) , the contour integration over u>\ can easily be per- 
formed. Taking the imaginary part leads to 

5&E ( _l(w > 0) = -U 2 / dep(e)377°(w-e) . (41) 
Jo 

The imaginary part of the second-order self-energy van- 
ishes for \w\ > 3W/2. The Hilbert transformation provides 
the real part as 

RI^M = -V / dC^i 2 ](C) N — . 

For practical calculations it is advisable to split the in- 
tegration routines into intervals [(r — l)W/2, rW/2] (r = 
1,2,3) in order to speed up the integrations and to mini- 
mize numerical errors. 

For the density of states of the Bethe lattice iJBJ, the 
result for the real and imaginary parts of the second-order 




self-energy are shown in Fig. ^ apart from the prefactor 
U 2 . As seen from the figure, the self-energy reproduces the 
Fermi-liquid relations (|14JI and H15|l for small frequencies. 
Explicitly, we find from our numerical integrations 

ZiU)- 1 = 1 + 1.307[1] +0{U 4 ) (43) 

and ^ 

7 (E/) = 3.242[1]Q^) +G(C/ 4 ), (44) 

where the number in brackets denotes the uncertainty in 
the last digit. 

3.2 Fourth order 

The twelve topologically different diagrams to fourth or- 
der can be grouped into four sets; see below. They were 
used earlier by Yamada and Yosida j22| in their study of 
the symmetric Anderson impurity model, and by Freericks 
and Jarrell |24| in their finite-temperature perturbation 
study of the Hubbard model. With the help of particle- 
hole symmetry it is not difficult to show that the diagrams 
of each set give the same contribution . We also find 

4 4) M = 3^( 4a )(c)+4 4b )(c)+4 4c )(c)+4 4d )H) . 

(45) 

We discuss the diagrams of the four sets and their contri- 
butions to the self-energy in the following. 



Fig. 2. Set A of three equivalent diagrams to the proper self- 
energy in fourth-order perturbation theory. 

are applied in the time domain, and particle-hole symme- 
try, G a {t) = — Gfj ( — t), is used appropriately. For the ring 
diagram the Feynman rules result in 



Ul-2 + Ull W 3 + Uli UJ,i + UJi 




(46) 

In order to evaluate this diagram, we define 

PH=[77»] 3 . (47) 
Analogously to the second-order calculation, we then find 

3l4 4a) (w > 0) = -t/ 4 / de /)(e)3P(w - e) (48) 
Jo 

with 

9fP( W ) = 3i7» [3(5ii7») 2 - (3i7») 2 ] . (49) 

The real part 5R7T°(w) is obtained via Hilbert transforma- 
tion of (uj); see appendix iBl 



3.2.1 Ring diagram 

The ring diagram in Fig. [21 and the particle- hole/particle- 
particle ladders give identical contributions at half band- 
filling. This is most easily seen when the Feynman rules 



3.2.2 Second-order diagram with second-order vertex 
correction 

The second-order diagram with second-order vertex cor- 
rection in Fig. [3 and the particle- hole/particle-particle 
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expressed in terms of help functions which can be tabu- 
lated and used in the remaining at most two-fold integra- 
tions over finite energy intervals. These help functions are 
listed in appendix IbI 

The ten remaining terms lead to 20 integrals which can 
be expressed using the bare density of states and the help 
functions. Appropriate changes in the energy integration 
variables allow us to regroup them into six terms, 




+ 



3JcS 4b) H = 




+ 




Fig. 3. Set B of three equivalent diagrams to the proper self- 
energy in fourth-order perturbation theory. 



ladders with crossed interaction lines give identical contri- 
butions at half band-filling. For the second-order diagram 
with second-order vertex correction, the Feynman rules 
result in 




UI3 U3 ■ 



Explicitly, 

4 4b) H = u 4 



(50) 



dw 3 



00 dwi 

-co 27ri J-oo 2?ri 

x G> -ui- uj 3 )G° a (w - co 3 )77> 3 ) ■ 

(51) 



The contour integration over u> 3 can be performed us- 
ing the spectral representation of the bare polarization 
bubble and of the local Green function ({3"5jl in (|5l|) . 
The remaining contour integral over u>x then results in 
24 terms, four of which are zero. When we restrict our- 
selves to lo > and focus on the imaginary part, another 
four terms vanish. A change in integration variables shows 
that six terms come in pairs. Thus, only ten different terms 
need to be evaluated [2SI ■ The calculation of the imaginary 
part reduces the seven-fold integrations over the bare den- 
sity of states to six-fold integrations. Fortunately, these 
integrations can be grouped so that the final result can be 



3^ 4b >( W >0)=7T[/ 4 ^/ 4 



(52) 



i=l 



with 
h 



„W pW 

dah(a) / dbh(b) 
/o Jo 

p(uj — a — b)l{b — oj)l(a — lj) , 

rW r W/2 



(53) 



h = - 



2 / dbh(b) / de x p{ei)h{uJ - e x ) 
Jo Jo 

[/(6-w)/(6-6i) + /(6 + w)/(6 + ei)] , (54) 



(55) 



h = 



nW pW/2 

2 / da h(a)p(ui — a) / deip(ei)/(ei + a) 
Jo Jo 
[H(ei +a-uj) + H(uj + e x )] , 

7r 2 / dah(a) / db h(b)p(ui — a)p(ui — b) 
Jo Jo 

\p(uj — a — 6) — p(a + b — lj)] , (56) 

pW/2 pW/2 

2 / deip(ei) / de 2 /9(e 2 )/i(a; - ei) 
Jo Jo 

fie, +e 2 -cj) [H(ex + e 2 ) + H(e 2 - lj)] , (57) 

W/2 pW/2 

deip(ei) / de 2 p(e2)p(oJ - e% - e 2 ) 
Jo 

2H(e 2 - uj)H{ei + e 2 ) + H(lj - e 2 )H(e 1 + e 2 ) 



Hie! - Lj)H{e 2 - lj) 



(58) 



where p(x) = p(x)0(x). For practical calculations, it is 
advisable to split the integration routines into intervals 
[(r - l)W/2, rW/2] (r = 1, . . . , 5) in order to speed up the 
integrations and to minimize numerical errors. We have 
checked our results against a numerical integration in the 
time domain which is easier to implement but much less 
accurate for the same computational effort. 



3.2.3 Second-order diagram with vertex correction in the 
polarization bubble 

The second-order diagram with vertex correction in the 
polarization bubble in Fig. 0] and the particle-hole/par- 
ticle-particle ladders with crossed interaction lines give 
identical contributions at half band-filling. For the second- 
order diagram with second-order vertex corrections in the 
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+ 




+ 




Fig. 4. Set C of three equivalent diagrams to the proper self- 
energy in fourth-order perturbation theory. 



polarization bubble, the Feynman rules result in 



0J4 — Cd2 



4 4c) (a 




(J — u± 



Obviously, we may split off the renormalized polarization 
bubble. Analogously to the second-order calculation, we 
immediately arrive at 



in pairs. Of the remaining 13 terms, five do not provide a 
finite imaginary part for x > 0. A re-grouping of integra- 
tion variables shows that only five different terms remain 
for x > 0, 

/>oo />oo 

3i7 v (x) = 23 / de 1(0 ( £l ) . . . / de 6 p(e 6 ) 
Jo Jo 



x + ei+ e 3 + e 4 + e 6 £1 + £3 + £4 + £5 
1 

x : 1 \ 

x + £1 - e 2 - ir?isgn(ei - e 2 ) 

1 1 

X - £2 - £3 - £4 - £5 + i?/ 3 £2 + £3 + £4 + £6 

1 



x + ei - e 2 - i?7isgn(ei - e 2 ) 

1 1 
; 

J a; - £1 - £ 2 + iryi .x - £ 2 - £3 - £4 - £5 + »/5 

1 



+2- 



£2 + £3 + £4 + £6 
1 1 



X - £1 - £2 + 1771 CC + £5 + £6 - 17/2 
1 

X 



£2 + £3 + £4 + £6 
I 1 



.E - £1 - £ 2 + l?7i X - £5 - £ 6 + 1772 
I 

X 

x - £ 2 - £3 - £4 - £5 + 1773 



(62) 



Taking the imaginary part simplifies this expression. First, 
the terms with the common factor sgn(£i — £2)<5(a; + £i — £2) 
cancel each other. Second, we are left with five- fold in- 
tegrations over finite intervals. Third, these expressions 
can be simplified further using the help functions of ap- 
pendix We find 



Q7T v (x > 0) = -2tt^ Ji 



(63) 



i=l 



with 



SE¥ c) (uj > 0) = -U 4 / dep(e)S*2Z'v(w - e) (60) J x = / d£ 2 / 

Jo Jo Jo 



with 

J-00 27Tl J -00 2lT1 
(3°(w4-w 2 )G^(« 4 -s-« a )IZ2(w a ) , (61) 

which is symmetric in x. 

Using the spectral representation of the bare polar- 
ization bubble H39fl and of the local Green function l|38|l 
in H6I(I , the contour integration over uj 2 results in six non- 
vanishing terms. The remaining contour integral over lu^ 
then gives 24 terms of which seven are zero, and four come 



W/2 nW 

d£ 2 / dy 
Jo 

h{y)p(e 2 )p(x - e 2 - y)f(x - e 2 )/(e 2 + y) , (64) 

J a = -2 / d£ 2 / dy 
Jo Jo 

h (y)p( e 2)p(x - e 2 )f(e 2 +y- x)f(e 2 + y) , (65) 

J 3 = -2 / d£ 2 / dy 
Jo Jo 

h(y)p( e 2)p(x - £ 2 - y)f(e 2 ~ x)f(e 2 + y) , (66) 

rW/2 rW/2 
Ji = 2 / d£ 2 / d£ 6 

Jo Jo 
p(e 2 )p(e 6 )p(x - e 2 )f{x + e e )H(e 2 + £e) , (67) 
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rW/2 pW/2 

J 5 = 2 / de 2 / de 5 

Jo Jo 
p(e2)p(e 5 )/o(x - e 2 )/(e 5 - x)H(t 2 + e 5 - x) , 

(68) 

J 6 = de 2 / de 5 

Jo Jo 

p(e2)p(e 5 )/(e 2 - x)f(e 5 - x)h{x - e 2 - e 5 ) , (69) 

pW/2 rW/2 

J 7 = -7r 2 / de 2 / de 5 
Jo Jo 
p(e2)p(e 5 )/o(x - e 2 )p(x - e 5 )h(x - e 2 - e 5 ) , (70) 

where p(x) = p{x)0{x). For practical calculations it is 
advisable to split the integration routines into intervals 
[(r - l)W/2, rW/2] (r = 1, . . . , 4) in order to speed up the 
integrations and to minimize numerical errors. We have 
checked our results against a numerical integration in the 
time domain. 

3.2.4 Second-order diagram with second-order self-energy 
insertion 




Fig. 5. Set D of three equivalent diagrams to the proper self- 
energy in fourth-order perturbation theory. 

The second-order diagrams with self-energy insertion 
in Fig. [5] give identical contributions in the paramagnet 
at half band-filling. These diagrams are not skeleton dia- 
grams so that momentum conservation cannot be ignored 
at the inner vertices. 



For the second-order diagram with self-energy inser- 
tion, the Feynman rules result in 



^ 4d) H = 




sL^-ui) 

( 71 ) 

Obviously, we may use the results for the bare polarization 
bubble and the second-order self-energy to simplify the 
expression to the form 

J-oo 27T1 

i^[G° CT ( e (k);c- Wl )] 2 . (72) 

k 

We must evaluate 

i£[G°( e (k);.x)] 2 = A + (x) +A_(x) , (73) 

k 

with 

rW/2 l l 

A±(x) = / dyp(y) — — — . (74) 

Jo x t y ± m xTy±m 

We may set 773 = 774 later so that A± (x) may be expressed 
in terms of p(0) and the derivative of the bare density of 
states, d(x) = (dp(x))/(dx). 

The spectral representations of A±(x) in (|74|l . II® {uS) 
in ^ 128(1 . and S^\lo) in l|131|l allow us to perform the 
contour integral over u)% in l|72|l with the result 

W/2 

3^ 4d) (w) = ttl7 4 J dxd{x)h(u - x)[S{x) - S{-x)] 


p3W/2 

+7TU i dbs(b)h(iu-b)l'(b) , (75) 

Jo 

where l'(b) denotes the derivative of the help function 1(b), 
see l(124|l . Again, it is advisable to split the integration 
routines into intervals [(r — l)W/2, rW/2] (r = 1, ... ,5) 
in order to speed up the integrations and to minimize nu- 
merical errors. 

The result for the real and imaginary parts of the 
fourth-order diagrams are shown in Fig. apart from the 
prefactor [/ 4 . As seen from the figure, the diagrams are 
equally important. At half band- filling there is generally 
no reason to include only special diagram classes, as done, 
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where the number in brackets denotes the uncertainty 
in the last digit. The prefactor of the fourth-order term 
in [^(J/)] -1 is smaller than reported in [Hj. There, the 
momentum distribution was calculated using Rayleigh- 
Schrodinger stationary perturbation theory. We suspect 
that stationary perturbation theory is flawed when dia- 
grams with self-energy insertions as in 1)71(1 occur. 

Fig. [7| shows that the fourth-order contribution to the 
self-energy becomes positive for < u> < 1.1W. For 

small interaction strengths this positive contribution is 
compensated by the overall negative second-order self- 
energy. The sum of both terms is no longer negative for 
all frequencies above U = 0.6AW, in contrast to the exact 
result J57J. This limits the applicability of fourth-order 
perturbation theory to moderate interactions strengths. 



Fig. 6. Imaginary part of the fourth-order self-energy dia- 
grams; units W = At = 1. 



3.3 Ground-state energy and average double 
occupancy 



e.g., in the RPA or in the ladder approximation. In partic- 
ular, for small oj the contribution of the fourth-order ring 
diagram is to a large extent canceled by the two diagrams 
with vertex corrections. 




Fig. 7. Real and imaginary part of the fourth-order self- 
energy; units W = At = 1. 



Fig.[7|shows the real and imaginary parts of the fourth- 
order self-energy, apart from the prefactor U . The self- 
energy reproduces the Fermi-liquid relations (|14|l and l|15|l 
for small frequencies. Explicitly, we find from our numer- 
ical integrations 



Z,/')- J =1 + l..:-!()7[:i](^ +0.739[1]Q^ 



and 



-(D=8.242[:i](g: 



-0{U 6 ) 
(76) 



-4.22[1]QQ +0(f/ 6 ), (77) 



The ground-state energy can also be expanded in terms 
of powers in (U/W) 2 , 



Eo 
L 



2W 



aW 



U 



-pw 



u 



3vr \W J 1 \W, 

We use this expansion in (|24() and obtain 



Mi 



2W 
~3T 



3a W 



■5(3W 



■0(U 6 ) . (78) 



-0(U 6 ) . (79) 



On the other hand, we may use (|23|l into which we insert 
the density of states from l(21|l which we can express in 
terms of the self-energy with the help of i|17[l and l|18|) . 
For a better numerical accuracy, we calculate 



3W/2 



b = 



+ Vu 
5W/2 



(80) 



-r(. 4 )(o,) + V^-^ 2) H] -4 



+ E?> M - V [w - E?> {u)-E [ V{u)Y-A 

(81) 

We fit the results to the form a = a Q (U/W) 2 +a 1 {U/W) A 1 
and b = boiU/W) 4 + &i([//IV~) 6 , and obtain for a = 
ao/(37r) and (3 = 6q/(57t) 



a = -0.08346[1] , 
= -0.0062[2] . 



(82) 
(83) 



The prefactor of the fourth-order term, (3, is an order 
of magnitude smaller than the second-order prefactor, a. 
This is very similar to the result for the single-impurity 
Anderson model addressed by Yamada and Yosida |23j . 
Therefore, we expect that the ground-state energy is rea- 
sonably well-described by fourth-order perturbation the- 
ory up to U < W. 
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Starting from 
we obtain the fourth-order result 

= 0.25 - 0.16692[2] (^j - 0.0248[8] (^j (85) 
+0(U 5 ) 

for the average double occupancy at half band-filling. The 
average double occupancy is positive semi-definite. This 
criterion excludes the applicability of (|85|l for interaction 
strengths U > 1.12 [1]W. 

An even more stringent condition can be drawn from 
a comparison with the ground-state energy of the Mott- 
Hubbard insulator at half band-filling ^3]. Up to third- 
order in l/U, 



Eo 
L 



U 
A 



W 4 



w 2 



32C7 512C/ 3 



+ 0{U^), (86) 



and 



(87) 



If we use the 1/ [/-expansion down to U > W, we find that 
the average double occupancy from the expansion in U 
and in l/U become equal at the crossing point U cmss — 
1.056W. Note that in [TSj we estimated U cA = 1.105[10]W 
for the opening of the gap. The fact that U CIOSS ~ U Ct i 
might indicate that the Mott-Hubbard transition is indeed 
continuous at U c ,i = U c without a discontinuity in the 
average double occupancy. 

Within perturbation theory, the ground-state energies 
from weak and strong coupling do not become equal. Fur- 
ther terms in the l/U expansion are also negative, so that 
eq. ((86(1 appears to give an upper bound to the ground- 
state energy. On the other hand, the perturbation expan- 
sion in U gives a maximum in Eq(JJ) at U = 1.12[1]W and 
thus provides a lower bound around U — W. Therefore, 
comparing energies from perturbation expansions will al- 
ways result in a lower ground-state energy for the 'metal' 
(expansion in U) than for the 'insulator' (expansion in 
l/U). In the above case we find the minimum energy dif- 
ference at U cross , AE (U CIOSS )/L = 0.018[1]W. 

Further observables will be discussed in Sects. 01 and |SJ 



4 Dynamical Mean-Field Theory (DMFT) 

In this section, we summarize the Dynamical Mean-Field 
Theory and two of its numerical implementations. The 
self-consistent solution of the single-impurity Anderson 
model is carried out by the Fixed-Energy Exact Diagonal- 
ization (FE-ED) JSj and the Dynamical Density-Matrix 



Renormalization Group (DDMRG) (14115) . We compare 
our results for the density of states from perturbation 
theory with those from FE-ED, DDMRG, the Numerical 
Renormalization Group (NRG) ^Ojj and Iterated Pertur- 
bation Theory (IPT). 



4.1 Single-Impurity Anderson Model 

In the limit of infinite dimensions [Hj and under the con- 
ditions of translational invariance and convergence of per- 
turbation theory in strong and weak coupling, lattice mod- 
els for correlated electrons can be mapped onto single- 
impurity models (2612713) , which must then be solved self- 
consistently. In general, these impurity models cannot be 
solved analytically. 

For an approximate numerical treatment various dif- 
ferent implementations are conceivable. One realization is 
the single-impurity Anderson model in the 'star geome- 
try', 



#SIAM = E zi$t-,$°\t- + U (d+d^ - - j ( d+d l - i 

t = U<7 ^ 
n s -l 



where Ve are real, positive hybridization matrix elements. 
The model describes the hybridization of an impurity site 
with an on-site Hubbard interaction to n s — 1 bath sites 
without interaction at energies ei < ti < . . . < i- In 
order to preserve particle-hole symmetry, we must choose 
e n ,-« = — ££, and Vn a -i = Ve for 1=1,..., n s — l. For even 
n s , this implies that there is an energy level at e„ s /2 = 0, 
i.e., at the impurity level. For this bath site, which is ab- 
sent for odd n s , we expect particularly strong hybridiza- 
tion with the impurity level, so that large odd-even effects 
can be expected. 

For a given set of parameters (eg, Vi) the model 1(88)1 
defines a many-body problem for which the one-particle 
Green function 



T 



SIAM 



(89) 



must be calculated numerically. Here, (. . .)siam implies 
the ground-state expectation value within the single-im- 
purity Anderson model. The various implementations of 
the DMFT differ in the choice of this 'impurity solver' (28). 
denoted, e.g., as DMFT(FE-ED) or DMFT(DDMRG). 

Ultimately, we will be interested in the limit n s — > oo 
where the hybridization function 



W (n.) (w)= £ 



V 2 

1 e 



lu — ei + ir)sgn(uj) 



(90) 



is required to approach the hybridization function of the 
continuous problem smoothly, 



H(u) = lim n ins) (to) . 



n s — >oo 



(91) 
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Correspondingly, the Green function is required to fulfill 
G a (uj)= lim Gfr-\u). (92) 

n s — >oo 

At self-consistency the Green function of the impurity 
problem describes the Hubbard model in infinite dimen- 
sions. As shown in the hybridization function must 
obey the simple relation 

H(u>) = G„(u) (93) 

on the Bethe lattice. This equation closes the self-con- 
sistency cycle for the continuous problem: bath energies 
and hybridizations must be chosen in such a way that the 
one-particle Green function and the hybridization function 
fulfill gEJ. 

4.2 Fixed-Energy Exact Diagonalization (FE-ED) 

In the FE-ED ^21 f° r the metallic state, we choose to 
resolve the density of states within the frequency inter- 
val \u>\ < W* /2. To this end, we partition the effective 
bandwidth W* into equidistant intervals of width 5W — 
W*/(n s - 1), 

I t = [-W*/2+(£-l)6W,-W*/2+£SW] , 1 < £< n s -l . 

(94) 

The energy levels are now fixed at the centers of Ig, 

e e = -W*/2 + (£- 1/2)SW , 1 < I < n a - 1 . (95) 

For our calculations we choose W* = 9t = 2.25W which 
is sufficient to resolve the main features of the density of 
states (quasi-particle peak, Hubbard bands) for weak to 
intermediate coupling strengths, U < W. In the FE-ED, 
we can be sure that our resolution increases systematically 
as a function of l/n s which cannot be guaranteed in the 
Caffarcl-Krauth implementation of the exact diagonaliza- 
tion |3I9I12I29| . 

For n s < 15, we determine the impurity Green function 
using the (dynamical) Lanczos technique. The imaginary 
part of the Green function displays Ul peaks depending 
on the number til ~ 100 states kept in the Lanczos di- 
agonalization. We apply a constant broadening of width 

5W to the individual peaks in \&r), t = 1, . . . , til. 

We then collect the weight into the intervals h and assign 
this weight to we — . Typically, the weight of the peaks 
at energies |w r | > W* /2 is very small. Thus we set 

/riL 
dc^3G^>K) 
-* r=l 

0(uj -u r + SW/2) - 0{uj -u r - SW/2) 

m ■ (96) 

In order to generate an educated input guess for Vi — 
y/w£, we use the results from fourth-order perturbation 
theory in (|96|l . In the FE-ED, the self-consistency cycle is 
stable, i.e., different input choices result in the same self- 
consistent solution for a given n s . This is not guaranteed 



in the Caffarel-Krauth exact-diagonalization scheme, as 
shown for the Mott-Hubbard insulator in ^21 ■ The reason 
for the stability of the FE-ED is that the energies et are 
kept fixed. The n s substantial peaks in QG^J 1 ^ (ui r ) carry 
sufficient information to determine the n s hybridization 
matrix elements Vi . However, these n s peaks are not suf- 
ficient to determine 2n s parameters (q, Ve) uniquely. 



T 




Fig. 8. Hybridization function 7v- ns \ui) for n s — 14 (cir- 
cles) and n a — 15 (squares) in comparison with the density 
of states from fourth-order perturbation theory at U = 0.4VF 
(solid line); units W = 4t = 1. 

Fig. |H1 shows the converged hybridization function in 
FE-ED for n s = 14 and n s = 15 for U = OAW = 1.6*, 
in comparison with the results for the density of states 
from fourth-order perturbation theory. It is seen that the 
agreement is very good. 



T 




CO 



Fig. 9. Hybridization function TiS n3 '(u)) for n s — 14 (cir- 
cles) and n a = 15 (squares) in comparison with the density 
of states from fourth-order perturbation theory at U = 0.6VF 
(solid line); units W = At = 1. 
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As seen in Fig. El deviations are noticeable for U = 
0.6W = 2At. Unfortunately, the resolution SW = 9t/U » 
0.1614 7 is still rather limited, and finite-size effects are 
appreciable around u> = 0.5W. Therefore, it is unclear 
whether the differences in D a {u>) between FE-ED and per- 
turbation theory around ui — 0.5W are due to finite-size 
effects of the FE-ED, or due to missing sixth-order terms 
in the perturbation expansion. The DDMRG which is able 
to treat much bigger systems allows us to resolve this is- 
sue. 



4.3 Dynamical Density-Matrix Renormalization Group 
(DDMRG) 

The Density-Matrix Renormalization Group (DMRG) is 
a numerical method in which the energy functional 



E{\*)) = (#\H\1>) 



(97) 



is minimized variationally in the subspace of normalized 
states (&\&) = 1. This is done by carrying out a numerical 
diagonalization on a system of finite size and a renormal- 
ization on a part of this system. At each step of the renor- 
malization, m density-matrix eigenstates of the subsystem 
are kept to build a variational basis of dimension 0(m 2 ); 
for a review, see |3U| . The optimal state for a given m, 
|^)opt; provides a variational bound for the ground-state 
energy, £ , Q ar = pt{&\H\\P) pt- For one-dimensional lattice 
systems, this method provides a highly accurate estimate 
for ground-state properties of hundreds of interacting elec- 
trons. The accuracy of the energy, E™ — Eq + 0(e 2 ), 
is better than that of the variational ground state itself, 
|^)o P t = \&o) + 0(e). Typically, the error e 2 scales as the 
weight P m of the discarded density-matrix eigenstates. 

In DDMRG ^1], this concept is generalized to the min- 
imization of a frequency-dependent functional. For the lo- 
cal Green function and uj > 0, it reads 



W v (\¥y,w) = (¥\(Eo+u-H) + r, 2 \&) 
+77<<f r o|c i , (T |>f) +r)(&\c+ \& ) 



(98) 



Here the frequency u> is fixed for an individual DDMRG 
run. The optimal state |^)opt is the imaginary part of 
the so-called correction vector. The optimal functional be- 



E +lu - H 

2 



^ (E + lu - E n y + v 2 



l#o> 
(99) 



where \W n ) and E n are the exact eigenstates and ener- 
gies of the Hamiltonian. The spectral representation (25} 
shows that, in the thermodynamic limit and up to a fac- 
tor of —ttt], the DDMRG provides the exact local density 
of states at frequency to, convolved with a Lorentzian of 
width 77, 

D2( W ) = -J- W°*(uj). (100) 

7T?7 



The main advantage of this variational approach is that, as 
in the ground-state energy calculations, the optimal value 
of the functional (i.e., the density of states) is obtained 
with an accuracy of the order of e 2 if the optimal state 
I \P) opt is calculated with an accuracy e. 

As the DMRG method is most accurate for systems 
with a quasi one-dimensional structure, we perform calcu- 
lations of the single-impurity Anderson model (|88|l in its 
equivalent linear-chain form |31j 



1 



-ffsiAM = U [ djd^ - - I I djd l - - 



a 

n s -2 



The DDMRG provides the exact local density of states 
for a finite chain with n s sites. To obtain the spectrum 
of an infinite chain, the broadening 77 must be scaled as a 
function of the system size ^3] . If 77 is chosen too small, the 
DDMRG density of states displays finite-size peaks. If 77 is 
chosen too large, relevant information is smeared out. As 
an empirical fact, 77 « W* /n s should be chosen, i.e., the 
resolution scales as the inverse system size, as found for 
one-dimensional lattice models. The DDMRG technique 
for the single-impurity Anderson model will be described 
in more detail elsewhere |15j . 

In the DDMRG extension of the FE-ED, the limita- 
tions of the Lanczos technique are overcome in two ways. 
First, we may use bigger systems. The DDMRG can han- 
dle n s = 64 sites on a workstation (m = 300 states kept 
in the renormalization) with CPU time as the limiting 
factor (48 hrs on a 500 MHz DEC- alpha workstation). In 
contrast, the exact diagonalization studies are seriously 
limited by memory constraints (3 GByte of memory for 
n s = 15). 

Second, the DDMRG provides the density of states 
at selected frequencies Wj. Typically, we choose them to 
resolve the effective bandwidth W* equidistantly, — 
LOi = Suj w 77 w 5W. We then 'deconvolve' the DDMRG 
data by inverting the Lorentz transformation 



5u) 



V 



n r] 2 + {uji - uj 3 )- 



(102) 



The procedure can be repeated for different choices of the 
equidistant frequencies u>i to get more values of D a (u>j). 
In this way, the DDMRG provides a set of values D a (cjj) 
for the density of states which, for the same n s , is at least 
as dense as the set from the Lanczos diagonalization. In 
practice, we use from one to four different sets of frequen- 
cies, corresponding to a frequency resolution from 77 to 
77/4. Naturally, structures with an intrinsic width of less 
than 77 cannot be resolved with this procedure even if we 
use different sets of frequencies. The main advantage of 
this transformation, however, is that no extrapolation or 
scaling analysis of these values D a (ujj) is necessary be- 
cause they converge very quickly to the n s — > 00 limit. 
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Therefore, with DDMRG we obtain a more accurate dis- 
crete representation of the density of states for a given n s 
than within the Lanczos FE-ED. 




-0.5 0.5 

CO 

Fig. 10. Density of states for U = from DDMRG on n a = 64 
sites, compared with the bare density of states po(u>) © (full 
line); units W = it = 1. Squares: data from DDMRG for 
5uj = rj = O.lt. Circles: DDMRG result after deconvolution. 

An example for the non-interacting case is shown in 
figure where we have chosen n s — 64, W* = 5t, 
6W = W*/63, and 8w = r\ = O.lt. This is a relevant test of 
accuracy because the non-interacting single-impurity An- 
derson model poses a non-trivial problem to DDMRG |15| . 
As seen in figure there is an excellent agreement be- 
tween the 'deconvolved' numerical data for n s = 64 and 
the exact result for n s — > oo. 



t 1 1 1 1 1 1 1 1 1 r 




Fig. 11. Density of states for U = OAW from DDMRG on 
n s = 32 sites after deconvolution (circles), compared with 
fourth-order perturbation theory (full line); units W — At = 1. 

The density of states from DDMRG and fourth-order 
perturbation theory also agree very nicely for U = OAW, 



i 1 1 1 1 1 1 1 1 1 r 




Fig. 12. Density of states for U = 0.6W from DDMRG on 
n s = 32 sites after deconvolution (circles), compared with 
fourth-order perturbation theory (full line); units W = At = 1. 

as is seen in figure ITTI This result confirms the reliability 
of both the DDMRG and our perturbation theory at weak 
coupling. 

The agreement is not perfect for U = 0.6W, as seen 
in figure El The DDMRG is applicable at all interaction 
strengths, it is limited only by finite-size effects (resolution 
1] cx l/n s ). Therefore, the deviations must be attributed to 
the sixth-order contributions which are missing in our per- 
turbation expansion. The magnitude of the discrepancies 
around u = 0.7W is indeed consistent with corrections of 
the order 0.6 2 = 36% relative to the fourth-order term, 
which dominates the second-order term around these fre- 
quencies. Therefore, we are confident that the DDMRG 
provides reliable results for all interaction strengths within 
its resolution limitations. 



4.4 Numerical Renormalization Group (NRG) 

The Numerical Renormalization Group (NRG) is a tech- 
nique which aims to resolve accurately the self-energy and 
the density of states near ui = 0. To this end, an exact di- 
agonalization of the single-impurity model is performed 
at the n-th step of the renormalization procedure for a 
frequency-interval 7™ around u) = 0. The width of the 
frequency interval is cut in half at the next step of the 
renormalization procedure, starting at I® — W*. There- 
fore, features near u> — are resolved with exponentially 
increasing accuracy whereas higher frequencies are rep- 
resented by a few peaks which are broadened on a log- 
arithmic scale ^01- As a result of this logarithmic mesh 
for the bath energies eg, frequencies of the order of the 
Hubbard bands are resolved with a rather limited accu- 
racy. Note that weight from the broadened Hubbard bands 
enters the low-energy physics again through the iterative 
self-consistency procedure. 

Fig. H3lshows the comparison of the self-energies from 
NRG and fourth-order perturbation theory for U = 0.4W^. 
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Fig. 13. Real and imaginary parts of the self-energy from 
NRG for U — OAW (dashed line), compared with fourth-order 
perturbation theory (full line) I32|. Also shown is the second- 
order self-energy (dotted line); units W = At = 1. 



At intermediate to large frequencies, the NRG does not re- 
produce the fingerprints of the Hubbard bands as already 
seen in perturbation theory. The same quantitative differ- 
ences are found for the real part of the self-energy. 

As a consequence, the density of states in NRG does 
not show any signs of the Hubbard bands at U = OAW, in 
contrast to the exact result. Fig. 1141 shows the density of 
states from NRG as compared with perturbation theory 
for U = OAW and U = O.&W. Even at U = 0.6VF, where 
the Hubbard bands are clearly visible, the NRG displays 
only a small shoulder. We conclude that the NRG has 
problems at frequencies of the order of the Hubbard bands 
whose height and positions cannot be determined reliably. 



4.5 Iterated Perturbation Theory (IPT) 

The IPT approximation to the self-energy of the Z — > oo 
Bethe lattice is given by |34I35| 



£i pT M = u 2 [°° ^n{n)g a {u n) 



(103) 



Fig. 14. Density of states for U = OAW and U = 0.6W from 
NRG (dashed lines) E2, IPT (dotted lines) [33], and fourth- 
order perturbation theory (full lines); units W = At = 1. 



Here, G a {ui) is the host Green function 

1 



oj - G a {oS) 



and we have defined 



dwi 
27rT 



(104) 



(105) 



as the polarization bubble of the host Green functions. 

By construction, IPT reduces to second-order pertur- 
bation theory in the limit of weak coupling. However, it 
omits the fourth-order diagrams with vertex corrections 
shown in figures and 01 whereas it weighs the other di- 
agrams differently. Thus, it is uncontrolled at moderate 
interaction strengths. 

Fig. ^] shows the density of states from IPT for U — 
0.4I-F and U = 0.6W in comparison with fourth-order per- 
turbation theory. It is seen that the IPT does surprisingly 
well at weak coupling considering that it is exact only to 
0(U 2 ). In this limit and for half band- filling, it appears 
to be superior to the results from NRG. However, it does 
not reproduce the height of the Hubbard bands correctly. 



Florian Gebhard et al.: Fourth-Order Perturbation Theory for the Half-Filled Hubbard Model in Infinite Dimensions 15 



5 Random Dispersion Approximation 

In this section, we present numerical results from the Ran- 
dom Dispersion Approximation (RDA). It becomes exact 
for lattice electrons in high dimensions. In contrast to the 
FE-ED, DDMRG, or NRG, it is not based on the DMFT 
self-consistency equations. Moreover, the RDA does not 
require the convergence of the perturbation expansion. 
Therefore, it provides an independent check of the validity 
of perturbation theory and the DMFT approach. 



5.1 Method 

In the Random Dispersion Approximation, the dispersion 
relation e(k) in the kinetic energy is replaced by a random 
quantity e RDA (k) where the bare density of states acts as 



the probability distribution 



RDA 



(k)) 



(106) 



and all correlation functions factorize according to 
etc. This is the characteristic property of the dispersion 
relation in infinite dimensions |2I6| , so that the RDA with 
the semi-elliptic density of states (01 becomes exact for 
the Bethe lattice with infinite coordination number. 

In order to put this idea into practice, we choose a 
one-dimensional lattice of L sites in momentum space 



kt = 



2tt 

T 



L + l 



(107) 



and determine the dispersion relation e(k) as the solution 
of the implicit equation 

fc/2 = (2e(k)/W) [l-(2e(fc)/PF) 2 ] 1/2 + arcsin(2e(fc)/W r ) . 

(108) 

This choice guarantees p(e) = pv(t) m the thermodynamic 
limit. 

Next, we choose a permutation Q a for each spin di- 
rection a which permutes the sequence {1,...,L} into 
{Q<t[1], • ■ ■ i Qa [-£]}■ This defines a realization of the RDA 
dispersion, Q = [Qj, QJ. The numerical task is then the 
Lanczos diagonalization of the Hamiltonian 



# S = EE £ (WC 



a c k t ,a 



UD 



In this way we obtain the momentum distribution 
n Q (e;U)= 1 -Y: 



i{k t )-- 



(109) 



(110) 



where (. . .) denotes the ground-state expectation value for 
the realization Q. 

As a next step, we obtain all physical quantities for 
fixed system size L by averaging over Nq realizations Q. 
Typically, we choose at least Nq = 100 for 6 < L < 14, 
and Nq = 50 for L = 16. For the physical quantities, 



we obtain gaussian-shaped distributions for which we can 
determine the average values, e.g., 



n(e;U) 



(111) 



with accuracy 0{1/Nq). 

In order to improve the quality of our distributions 
slightly, we impose a filter on our randomly chosen permu- 
tations. For a truly random dispersion, L|i RDA (£)| 2 = e 2 
is independent of £ [2J- Therefore, we discard those real- 
izations for which 



L-l 

E 

e=i 



[\t Q °(l)f 



> d L 



(112) 



with cLl ~ 0.2 for L < 16. In this way, we admit about 
every second of the randomly chosen configurations. Note 
that there are of the order of (LI) 2 different realizations 
so that our filter does not introduce any unwanted bias. 



5.2 Momentum distribution 
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Fig. 15. Momentum distribution in RDA for U = W/n ~ 
0.32W, U = 2W/n » 0.64W, and U = 5W7(2tt) w 0.80W for 
even system sizes L < 16. The lines are the result from fourth- 
order perturbation theory. Note the weak elne behavior near 
the discontinuity. 



Fig. 1151 shows the RDA momentum distribution in 
comparison with results from perturbation theory. As seen 
from the figure, finite-size effects are rather small in the 
metallic phase. All values for the momentum distribution 
appear to fall onto almost the same curve. For small in- 
teraction strengths, the RDA data lie essentially on top 
of the perturbative results. Deviations for larger interac- 
tions can be attributed to missing higher-order corrections 
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to our fourth-order result. Nevertheless, it is quite aston- 
ishing that fourth-order perturbation theory provides a 
sensible description of the momentum distribution for in- 
teraction strengths as large as U w 0.8VF. 



1 


I 1 1 1 


1 1 1 1 














— PT4 










IPT 










— ■ NRG 




0.5 











1 




-0.4 -0.2 ( 


) 0.2 


0.4 



e 

Fig. 16. Momentum distribution for U = 0.6W in IPT (dashed 
line) 33;, NRG (dotted line) [3*^1 , and fourth-order perturba- 
tion theory (full line). The deviations of NRG near the jump 
result from a small but finite positive imaginary part of the 
NRG self-energy around uj = 0. 

In figure 111)1 we compare our results for the momentum 
distribution at U — Q.QW from fourth-order perturbation 
theory with those from IPT and NRG. We see that both 
approaches describe the momentum distribution very well 
despite their shortcomings for the self-energy at interme- 
diate frequencies. Apparently, for small to intermediate 
coupling strengths, the overall momentum distribution is 
not a very sensitive test for the quality of an approxima- 
tion to the self-energy at intermediate frequencies. 

A more sensitive quantity is the discontinuity of the 
momentum distribution at e = 0, see (1301) . In figure El we 
show the quasi-particle weight as a function of U /W as 
obtained from IPT, NRG, RDA, and fourth-order pertur- 
bation theory. For the latter quantity, we invert eq. H76J) 
which gives 

Z(U) = 1 - 1.307[1] +O.969[2]Q0 + 0(U e ) . 

In the region where fourth-order perturbation is reliable, 
U < 0.6W, the results agree with those from NRG and 
RDA within their numerical accuracy. Therefore, fourth- 
order perturbation theory cannot discriminate in favor 
of either of the two approaches which, however, support 
different scenarios for the Mott-Hubbard metal-insulator 
transition. 

The quasi-particle weight from Iterated Perturbation 
Theory closely follows the fourth-order result to a point 
where (|113l) definitely overestimates Z(U), e.g., for U — 




0.5 1 1.5 2 

Fig. 17. Quasi-particle weight as a function of U/W in IPT 
(dotted line) [331, NRG ( nlled circles) [3J], RDA (squares), and 
fourth-order perturbation theory (full line). 

0.7W. From this we conclude that IPT generally overesti- 
mates the stability of the metal. 



6 Conclusions 

In this work we have calculated the one-particle self-en- 
ergy for the half-filled Hubbard model in the limit of infi- 
nite dimensions to fourth order in the interaction strength. 
We have reduced the four non-equivalent fourth-order dia- 
grams to sets of two-dimensional integrals over real, tabu- 
lated functions over finite intervals. From the self-energy, 
we have derived the density of states, the ground-state 
energy, the mean double occupancy, the momentum dis- 
tribution, and the quasi-particle weight, and have com- 
pared it to various analytical and numerical approaches 
to the Hubbard model in infinite dimensions for the case 
of a Bethe lattice with infinite coordination number, i.e., 
a semi-circular bare density of states. 

Our results show that, at half filling, it is not permitted 
to select only subclasses of diagrams because all diagrams 
contribute equally. The imaginary part of the self-energy 
to fourth order becomes positive at some frequency for 
U = Q.64W, which limits the applicability of our fourth- 
order results to moderate interaction strengths. For the 
momentum distribution and, in particular, the ground- 
state energy and the mean double occupancy, the pertur- 
bation expansion appears to be better behaved, i.e., the 
results remain sensible up to (/ i=s W. A similar obser- 
vation had been made for the single-impurity Anderson 
model by Yamada and Yosida. 
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The Dynamical Mean-Field Theory (DMFT), which 
becomes exact in infinite dimensions, requires the self- 
consistent solution of a single-impurity Anderson model. 
Our results indicate that the Dynamical Density-Matrix 
Renormalization Group (DDMRG) is an excellent 'impu- 
rity solver'. Where Exact Diagonalization is hampered by 
finite-size effects (n s < 15), the DDMRG treats much big- 
ger systems, up to n s = 64, and a new scheme for the 
deconvolution of the data further improves the frequency 
resolution. In this way, we obtain a very good agreement 
between the results from DDMRG and fourth-order per- 
turbation theory where the latter is applicable. 

The agreement is found to be less satisfactory with the 
results from Numerical Renormalization Group (NRG) at 
intermediate frequencies. The Hubbard bands which be- 
come discernible around U = OAW are not resolved very 
well within the NRG scheme. It would be interesting to 
see whether a deconvolution of the NRG data could im- 
prove the agreement with fourth-order perturbation the- 
ory at the energy scale of the Hubbard bands. Iterated 
Perturbation Theory (IPT) works surprisingly well for 
U < 0.6 for all quantities tested. However, for larger in- 
teraction strengths, it seriously overestimates the quasi- 
particle weight and thus the stability of the metallic state. 

The Random Dispersion Approximation (RDA) , which 
also becomes exact in infinite dimensions, provides the mo- 
mentum distribution and the quasi-particle weight. The 
results agree with those from fourth-order perturbation 
theory within the limits of its applicability. Unfortunately, 
for U < W, the momentum distribution and the quasi- 
particle weight turn out not be very sensitive quantities, 
i.e., all methods equally well reproduce the results from 
fourth-order perturbation. However, the RDA and NRG 
support two different scenarios for the Mott-Hubbard me- 
tal-insulator transition. The differences between NRG and 
RDA in the quasi-particle-weight become sizable only for 
U > 0.8W, a region which, unfortunately, cannot be ac- 
cessed with fourth-order perturbation theory. 

We hope that the DDMRG approach will enable us to 
investigate the region U > O.QW in more detail. Work in 
this direction is in progress. 
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A Momentum distribution near the 
discontinuity 

Using the Luttinger relations H14[) and l|15|l for a Fermi 
liquid, we evaluate J2B for n a (e > 0) = (1 - Z(U))/2 + 



Sn(e). Then, 
5n a (e) = 



did 



70; 



(114) 



7T (lj/Z - e) 2 + 7 2 w 4 

apart from regular terms which come from the incoherent 
background in l|25[) . Naturally, we are interested in the 
behavior close to the Fermi energy, i.e., for energies which 
fulfill 

< Z(U)e < uj c . (115) 

Therefore, we require that 

Z(U)qcuj c (116) 

in order to observe a finite e-interval in which the formulae 
below are applicable. Furthermore, we choose e > so that 
we are above the jump in the distribution. 

As can be seen by an explicit calculation, the dominant 
contribution for small e comes from the frequencies u) = 
Z(U)e. Then, we may ignore the term proportional to lu 4 
in the denominator and we find 



5n a {Z(U)e « w c ) = 



jZ(U) 2 lu c 



dx : 



This approximation requires that 

1 



! (x - Z(U)e/iu c ) 2 ' 
(117) 



Z(C/) 2 7 ■ 

This condition should correspond to l|115|l so that 
1 



oc 0J C 



(118) 



(119) 



Z(Uh 

7 oc M~ 2 . (120) 

The scaling laws (|116[) and i|120|) appear to be natural, 
given the Fermi liquid equations l|14f> and l|15fl . To leading 
order in elne, eq. (|117fl results in (|31|l . 

If the above scaling laws (|116fl and l|120[) apply, there 
should be a finite region around e ~ in which the e In e 
dependence can be seen. However, its overall intensity 
scales down with an extra factor Z(U). 

B Definition of help functions 

Here we define the help functions which are necessary for 
the evaluations of the diagrams to fourth order. Explicit 
results are given for the semicircular density of states © 
in units of W = 4i = 1 . 

We start with two functions which characterize the 
non-interacting Green function. 

f W/2 P(e) 



(121) 



ituj - 



1 - v 7 ! - 4w 2 In 



2w 



1 + vT - 4w 2 



t/(M > o) = ™ ~ 1 ~ sgnMv^ 



1 arccos 



(122) 
1 

(123) 
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8lu 1-0(\uj\ -1/2) Jl 



(4w 2 



(124) 



(125) 



Next, we express the imaginary part of the bare polariza- 
tion bubble in the form (0 < a < W) 

pW/2 pW/2 

h{a) = / dei / de 2 5 (a - ei - e 2 ) p(ei) p(e 2 ) 
Jo Jo 

(126) 

W-a 



w r 

+6{a- —) J dyp 



a+y\ a-y 



and 



r w , 1 

n» = - dah(a)( — 

Jo \uj-a + i 



(128) 



ir] cj + a — it] 
Moreover, we need the Hilbert transform of h(a) 

H(x)=V [ da^- , (129) 
Jo x + a 

so that 37I°(w > 0) = irh(uo) and 9W7°M = H(u) + 
H{-u)._ 

The imaginary part of the second-order self-energy re- 
quires the function 

nW/2 rW/2 pW/2 

s(b) = / dei / de 2 / de 3 
Jo Jo Jo 

p(ei)p(e 2 )p(e 3 )S(b - a - e 2 - e 3 ) , (130) 

which leads to 

r 3W/2 f , 1 x 

= U 2 / dbs(b) L —— + ; 

a ' Jo yJ \uj + b-ir] uj-b + irij 

(131) 

compare (|41|> . Moreover, we need the Hilbert transform of 

8(b) 

f 3W/2 s(b) 
S(x)=V d6^-, (132) 

Jo x + o 

so that 9Zi 2) (w > 0) = -ttU 2 s(uj) and V.S {2) {uj) = 
U 2 [S{tu)- S(-w)]. 

Finally, we introduce the derivative of the bare density 
of states as 



d(x) 



dx 

16 x 



tt VT - Ax 2 



for |x| < 1/2 . (133) 
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